--- title: Outlier Detection and Missing Imputation date: 2021-08-18 11:11:00 ---
최근 많은 분야에서 자동화가 이루어지고 있고, 자연현상의 관측 역시 관측장비의 발달로 자동화가 상당부분 이루어졌다. 생,화학분야는 상대적으로 미흡하지만 조석, 파랑, 수온 등 과 같은 해양 물리자료의 관측은 자동관측 시스템이 대부분이다.
자동화된 관측장비에서 생산되는 시계열 자료는 기존에 사람이 직접 관측하던 것보다 고해상도이기 때문에 기존보다 활용도가 높을 수 있으나, 기술적인 한계로 인한 결측과 이상값, 자연현상의 복잡함을 그대로 기록하는 특성으로 인해 원시자료를 그대로 사용할 때 분석 결과가 왜곡될 위험이 있다.
Acuna and Rodriguez(2004)는 분류모델의 경우, 1%미만: 심각한 영향없음, 1~5%: 관리 가능, 5~15%: 정교한 방법론 필요, 15%이상: 예측 모델의 심각한 성능저하 로 결측치 비율에 따른 결과 왜곡 현상이 있음을 언급하고 있으며, 강필성(2012)은 결측치를 제외하는 방법보다 간단한 선형 대치 기법을 통해서 자료를 보완하는 것이 분류모델의 성능저하가 적다고 언급했다.
이상자료의 경우도 결측과 마찬가지로 자료의 통계량을 계산하거나 모델을 구성할 때 왜곡을 발생시킬 수 있으며(조홍연 등 2016) 자료에 실제 존재하는 극값과의 구분이 어려워 실제 자료분석에서의 이상자료 탐지 및 제거기법의 활용은 미흡하고 경험에 의존하여 제거하기도 한다.
따라서 본 연구에서는 기상청에서 제공하는 해양기상관측부이의 수온자료를 이용하여 결측치의 대치와 이상자료 탐지 및 대치 기법에 대해서 논의한다.
wt <- kma_buoy_subset$Water_Temperature
## summary data
print(summary_wt <- wt %>% describe)
## .
##
## 1 Variables 87672 Observations
## ---------------------------------------------------------------------------
## Water_Temperature
## n missing distinct Info Mean Gmd .05 .10
## 83891 3781 253 1 17 5.438 10.8 11.6
## .25 .50 .75 .90 .95
## 13.1 15.6 20.9 24.0 25.6
##
## lowest : 0.0 0.1 0.2 5.8 5.9, highest: 30.7 30.9 31.3 31.4 31.6
## ---------------------------------------------------------------------------
## missing ratio = 4.5%
print(paste0("missing ratio = ", 100 * round(as.numeric(summary_wt$Water_Temperature$counts[2])/as.numeric(summary_wt$Water_Temperature$counts[1]), 4), "%"))
## [1] "missing ratio = 4.51%"
autoplot(kma_buoy_subset$Water_Temperature)
가장 먼저 기계적 오류로 인한 이상치로 판단되는 값부터 제거(과정 생략) 일반적으로 특정 값을 선택적으로 제거하거나(ex -999),rosnerTest같은 분포기반 이상치 탐지 기법으로 제거.
rosner test로 제거하는 과정이 계산시간 측면에서 효율적이지 못하여(코드 수행시간이 오래걸린다.) 본 과정에서는 해당 코드를 실행하지 않는다.
Outlier는 1.5IQR 범위를 벗어나는 mild outlier와 3IQR 범위를 벗어나는 Extreme Outlier로 구분할 수 있음.
여기서는 사전에 정의한 IQR_filter()함수를 이용하여 3IQR(default) 범위를 넘어서는 Extreme Outlier만 탐지하도록 설정.
IQR_filter
## function (x, type = "extreme")
## {
## if (type == "extreme") {
## IQR_range <- 3
## }
## else if (type == "mild") {
## IQR_range <- 1.5
## }
## else {
## stop("Error, correct type should be required")
## }
## q <- quantile(x, c(0.25, 0.5, 0.75), na.rm = T)
## lower <- q[1] - IQR_range * IQR(x, na.rm = T)
## upper <- q[3] + IQR_range * IQR(x, na.rm = T)
## return(ifelse(x < lower | x > upper, NA, x))
## }
product1 <- kma_buoy_subset$Water_Temperature %>%
apply(2, IQR_filter) %>% as.xts
p1 <- autoplot(product1)
p1
원 자료에 그대로 IQR filter를 적용해도 2009년 근처의 이상치 제거가 되지 않음.
Rosner Test는 코드 수행시간이 오래걸리기 때문에 현재 자료는IQR filter()함수 사용이 적합하다고 판단됨.
일괄 처리 문제가 있으므로 Shifting method를 사용.
imputeTS package이상치가 있는 구간에서 결측이 함께 발생할 경우 결측값 추정이 부적절하며,
따라서, 이상치에 강건한
LOWESS(Locally WEighted Polynomial Scatterplot Smoothing)기법을 이용하여 이상치를 탐지 및 제거한다.
또한, Short term missing에 대한 대치 성능은 크게 문제가 없어보이나, 다항식을 이용한 적합 시, 장기결측 구간에서 결측치 대치 성능이 크게 떨어지는 현상이 발생함.
LOWESS 적합과정에서 사용된 Bandwidth 는 전체 자료에 대해서는 적절한 Smoothing 으로 판단되지만, 부분별로 확대해서 살펴보게 되면 연별 변화양상보다 짧은 주기의 변동에 대해서는 Over Smoothing 이 된 것을 확인할 수 있다. 따라서 대부분 주기성분을 갖는 해양 관측자료를 Smoothing 할 때에는 가장 짧은 주기의 변동성분을 대표할 수 있는 Bandwidth 설정이 필요하며 특히 장기관측자료를 한번에 처리할 경우 단주기의 정보손실을 고려해야 한다.
아래 코드를 통해 잔차의 특성을 살펴본다.
## residuals1
##
## 1 Variables 87672 Observations
## ---------------------------------------------------------------------------
## Water_Temperature
## n missing distinct Info Mean Gmd .05
## 83891 3781 83891 1 -0.005396 1.292 -1.81137
## .10 .25 .50 .75 .90 .95
## -1.38475 -0.71024 -0.05238 0.66364 1.45776 2.02038
##
## lowest : -24.969102 -19.940401 -19.136007 -19.029932 -19.027409
## highest: 5.961000 6.261774 6.317315 6.418865 6.935617
## ---------------------------------------------------------------------------
Step 3-a 에서 imputeTS package의 na.interpolation() 함수를 사용하여 Raw data의 결측구간의 보간을 시도함. 짧은 구간의 보간은 비교적 잘 되었으나 장기결측 구간의 보간은 값을 직접 사용하기에는 문제 발생.
따라서 Step 3-b 에서 수행한 Locally Polynomial Smoothing 으로 계산된 값을 이용하여 다음 단계인 Step 4 에서 Step 1의 IQR filter Step 3-a 의 Interpolation과 EMD, SSA 기법을 적용해본다.
Step 3 에서 계산된 Smoothing Data 에 IQR filter를 다시 적용하며 이는 Masked Outlier 를 탐지해 내기 위함이다. 여기서는 정보손실을 최소화 하기 위하여 Extreme Outlier 탐지 조건인 3*IQR 범위를 적용하였으나, Over Smoothing 된 값으로 인해 Raw Data 의 정상치로 추정되는 값이 이상자료로 판단되어 제거되는 것을 방지하기 위해 더 완화된 3*IQR 이상의 조건을 고려해 볼 수 있다. 이는 자료를 Smoothing 할 때 발생하는 과적합 문제를 해결하기 위해 임의로 Bandwidth를 넓게 조정한 것에 대한 조치이다.
Raw Data 와 Smoothing Data의 잔차를 이용하여 이상자료를 탐지, 제거하고 난 후의 plot.
연별 변화패턴에 Masking 되어 있던 이상자료 중 극단적인 값들이 제거된 것을 볼 수 있으며 이로 인해 다시 계산한 Smoothing line 이 달라짐.
## [1] 6740
## [1] 254
## [1] 28
1.5 IQR, 3 IQR, 5 IQR일 경우, 각각 6740, 254, 28개의 이상치가 탐지됨.
결측 구간을 제외한 전체 자료 83891개 중 0.3% 미만의 극히 일부자료이기 때문에 본 분석에서는 3 IQR 조건과 5 IQR 조건의 차이는 미미하다고 할 수 있고 따라서 3 IQR조건을 적용하여 이론적인 근거(Extream Outlier)를 확보.
단, 1.5 IQR범위에서 필터링 할 경우, 전체 자료의 약 8%가 이상자료로 감지되기 때문에 고도화된 대치(Imputation)방법 강구 필요.
Step 1~4 과정을 통하여 Outlier 제거와 Missing Data 의 대치는 동시에 수행되어야 한다는 것을 알 수 있음.
따라서 Step 5 에서는 Step 4 까지의 Outlier 제거과정을 거치며 얻은 Preprocessed Data 를 이용하여 그럴듯 하게 Missing Data Imputation 처리. 각 기법의 적용은 imputeTS package 의 na.kalman(), na.locf(), na.ma(), na.seadec() 함수들과 그 인자를 조정해가며 테스트해본다.
na.seadec() 함수는 자료의 seasonality를 계산하지 못하여 사용불가.
정량적인 Imputation 성능평가는 어렵지만 장기 결측 구간의 대치 결과
na.kalman()함수의 성능이 비교적 나아보이며,na.seadec()함수는 장기 결측구간의 존재로 인해 주기성분의 계산이 불가능하여 오류가 발생하는 것으로 보인다. 현재 이를 해결하기 위해 앞서 언급한na.kalman()함수로 계산된 자료를 전처리 자료에 대치시키고na.seadec()함수를 다시 적용할 계획이다.
Missing Data 관련 오류 발생으로 Rssa::gapfill() 함수를 비롯한 Missing imputation 연산 실행이 되지 않았다. 결측구간이 지나치게 길어 프로그램이 다루지 못해 발생하는 문제일 것으로 추정된다.
따라서 spectral.methods 패키지의 gapfillSSA()함수를 사용했다.
series.ex <- as.numeric(raw3)
indices.gaps <- is.na(series.ex)
## perform gap filling
data.filled <- gapfillSSA(series = series.ex, plot.results = FALSE, open.plot = FALSE)
## 2018-08-24 00:02:08 : gap filling: Performing data preparation...
## 2018-08-24 00:02:09 : gap filling: Performing SSA iteration...
## 2018-08-24 00:04:40 : gap filling: Finished 1 of 10 outer iterations.
## 2018-08-24 00:06:51 : gap filling: Finished 2 of 10 outer iterations.
## 2018-08-24 00:09:05 : gap filling: Finished 3 of 10 outer iterations.
## 2018-08-24 00:10:35 : gap filling: Finished 4 of 10 outer iterations.
## 2018-08-24 00:12:05 : gap filling: Finished 5 of 10 outer iterations.
## 2018-08-24 00:13:36 : gap filling: Finished 6 of 10 outer iterations.
## 2018-08-24 00:15:07 : gap filling: Finished 7 of 10 outer iterations.
## 2018-08-24 00:16:38 : gap filling: Finished 8 of 10 outer iterations.
## 2018-08-24 00:18:10 : gap filling: Finished 9 of 10 outer iterations.
## 2018-08-24 00:19:42 : gap filling: Finished 10 of 10 outer iterations.
## 2018-08-24 00:19:42 : gap filling: Reconstructing timeseries...
## 2018-08-24 00:21:58 : Reconstruction: finished 1 of 7 outer iterations.
## 2018-08-24 00:24:17 : Reconstruction: finished 2 of 7 outer iterations.
## 2018-08-24 00:26:37 : Reconstruction: finished 3 of 7 outer iterations.
## 2018-08-24 00:28:09 : Reconstruction: finished 4 of 7 outer iterations.
## 2018-08-24 00:29:42 : Reconstruction: finished 5 of 7 outer iterations.
## 2018-08-24 00:31:14 : Reconstruction: finished 6 of 7 outer iterations.
## 2018-08-24 00:32:46 : Reconstruction: finished 7 of 7 outer iterations.
## 2018-08-24 00:32:46 : Reconstruction: Iteration process finished.
## plot series and filled series
impute_ssa <- xts(data.filled$filled.series, order.by=index(raw3))
impute_ssa[!is.na(raw3),] <- NA
binded_impute_ssa <- cbind(raw3, impute_ssa)
dygraph(binded_impute_ssa, main="raw & imputed data: SSA") %>%
dyAxis("y", label = "Temp (C)", valueRange = c(0, 40)) %>%
dySeries(label="raw") %>%
dySeries(label="SSA", strokeWidth = 1.5) %>%
dyOptions(colors = c("gray","brown")) %>%
dyRangeSelector()
SV <- read.csv("WT_obs.csv", header = T)[,-1]
jd <- julian( ymd(substr(colnames(SV), 2, 11)), origin = ymd(20150101))
p_dec2 <- function(x, jd, m){
SV <- as.numeric(x)
Tj <- 365.25/m
wj <- 2*pi/Tj
constant_term <- rep(1, length(jd))
cos_term <- t(apply(as.matrix(jd), 1, function(x) cos(wj*x)))
sin_term <- t(apply(as.matrix(jd), 1, function(x) sin(wj*x)))
colnames(cos_term) <- paste0("coef_", "cos", m)
colnames(sin_term) <- paste0("coef_", "sin", m)
PP <- cbind(constant_term, cos_term, sin_term)
CC <- solve(t(PP)%*%PP)%*%(t(PP)%*%SV)
eSV <- PP%*%CC
res <- list(CC, eSV)
return(res)
}
test_sig1 <- p_dec2(SV[1,], jd = jd, m = c(1,2))
test_sig2 <- p_dec2(SV[1,], jd = jd, m = c(1,2,3))
par(family = "serif", mar = c(5, 7, 5, 5))
plot(test_sig1[[2]], type = "l", lwd = 2, ylim = c(0, 30),
xlab = "Julian Days", ylab = "Amplitude",
main = "Combined Trigonometric Function",
las = 1, cex.main = 2, cex.lab = 2, cex.axis = 1.7)
lines(as.numeric(test_sig2[[2]]), lwd = 2, col = 4, lty = 2)
points(as.numeric(SV[1,]), pch = 16, cex = 0.5, col = 2)
RMSE1 <- sqrt(sum((test_sig1[[2]]-SV[1,])^2)/length(SV[1,]))
RMSE2 <- sqrt(sum((test_sig2[[2]]-SV[1,])^2)/length(SV[1,]))
AIC1 <- 2*nrow(test_sig1[[1]])-2*log(RMSE1)
AIC2 <- 2*nrow(test_sig2[[1]])-2*log(RMSE2)
c(k_5 = AIC1, k_7 = AIC2)
## k_5 k_7
## 9.597281 13.833785
삼각함수 조합 수를 최적화 하기 위한 기준으로
AIC(Akaike Inforamtion Criteria)를 사용.
AIC를cost function으로 사용하는 최적화 알고리즘을 결측자료 추정곡선 계산에 내장하여over fitting문제를 내부에서 해결할 수 있도록 설계.
RNN과 같은 신경망 구조에embedding할 수 있는지 검토 필요.
“Acuna, E., & Rodriguez, C. (2004). The treatment of missing values and its effect on classifier accuracy. In Classification, clustering, and data mining applications (pp. 639-647). Springer, Berlin, Heidelberg.”
“강필성. (2012). 분류 성능 향상을 위한 지역적 선형 재구축 기반 결측치 대치. 대한산업공학회지, 38(4), 276-284.”
“조홍연, 이기섭, & 안순모. (2016). 이상자료가 연안 환경자료의 통계 척도에 미치는 영향. Ocean & Polar Research, 38(2).”
“dplyr”, “ggplot2”, “Hmisc”, “zoo”, “tsoutliers”, “ggfortify”, “xts”, “EnvStats”, “KernSmooth”, “dygraphs”, “imputeTS”, “spectral.methods”
IQR(Inter Quartile Range)
locally polynomial regression
interpolation - linear, spline, stineman, moving average, kalman(arima model, structural model)
SSA(Singular Spectrum Analysis) based imputation